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Finding the global minimum of a multivariate function efficiently is a fundamental yet diffi- 
cult problem in many branches of theoretical physics and chemistry. However, we observe that 
there are many physical systems for which the extremizing equations have polynomial-like non- 
linearity. This allows the use of Algebraic Geometry techniques to solve these equations com- 
pletely. The global minimum can then straightforwardly be found by the second derivative test. 
As a warm-up example, here we study lattice Landau gauge for compact U(l) and propose two 
methods to solve the corresponding gauge-fixing equations. In a first step, we obtain all Gribov 
copies on one and two dimensional lattices. For simple 3x3 systems their number can already be 
of the order of thousands. We anticipate that the computational and numerical algebraic geometry 
methods employed have far-reaching implications beyond the simple but illustrating examples 
discussed here. 
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1. Introduction 

Finding the global minimum of a multivariate function is one of the most important tasks in 
statistical mechanics, condensed matter theory, lattice gauge theories and theoretical chemistry. 
Most of the usual methods to minimize a multivariate function are based on the Newton-Raphson 
method, where a start solution is guessed and then is refined by successive iterations. There are 
several efficient refinements of this minimization procedure, e.g., simulated annealing. However, 
they are often not successful as they can get trapped in a local minimum instead of the global one. 
This problem introduces an error of an unknown order. 

A typical example is the Landau gauge on the lattice: in the continuum the corresponding 
Landau gauge fixing condition d^A^ = is solved to fix the gauge, whereas on the lattice, usually 
the so-called standard lattice Landau gauge-fixing (LLG) functional (whose first derivatives with 
respect to the gauge parameters are the corresponding gauge-fixing conditions on the lattice) is 



numerically minimized [|9|, [lfj, |11], |12fl . We take the example of the LLG functional for compact 
U(l), which is given by 

f = £(i - cos(^ + e,- +A - (l.i) 

where the 0, £ (— 71, 7r] variables sit on the ith lattice site of a <i-dimensional lattice grid, jX is the 
unit vector in pL direction, and the 1)jU £ (—71, 7i] variables sit on the links between the rth and 
(/ + /i)th lattice sites. Here the 0-variables are related to the gauge potential in the corresponding 
continuum gauge theory and the 0-variables are the gauge transformations with respect to which 
F is minimized. The special case when all ^ are zero is called the trivial orbit. We will also 
frequently use random choices for the link angles 0, ijU which we then refer to as random orbits. The 
Hessian matrix of this functional is the Faddeev-Popov operator and its determinant is known as 
the Faddeev-Popov determinant JT7| ] on the lattice. 

In order to compare the results from lattice Landau gauge to those from functional methods in 
the continuum such as studies of Dyson-Schwinger equations [f7|], one should in principle find all 
stationary points of F under variation with respect to the gauge variables 0, [Q, [3| p, that is all 
solutions to the corresponding LLG equations on the lattice, 

dF 

M =0 - <12) 

The solutions to these equations are the lattice analogues of Gribov copies Ira. 1 The sum over all 
those Gribov copies weighted by the sign of the Faddeev-Popov determinant evaluated at each copy 
would then be independent of the gauge orbit by the Poincare-Hopf theorem and this could thus be 
used to define the measure of Landau gauge on the lattice analogous to standard Faddeev-Popov 
theory with Becchi-Rouet-Stora-Tyutin (BRST) symmetry. 

The reason why this procedure fails is that the signs of all Gribov copies exactly cancel one 
another, always yielding an exact zero for their sum which computes the vanishing Euler charac- 
teristic of the lattice gauge group manifold in the standard implementation of the Landau gauge. 
This zero sum is the origin for the famous Neuberger 0/0 problem of standard BRST formulations 
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on the lattice [O, 16], which is the statement that the expectation value of any gauge-invariant 



observable is of indefinite form 0/0 in such a formulation. 

The 0/0 problem can be overcome by stereographically projecting the lattice gauge group onto 
a manifold whose Euler characteristic is non-zero This provides a well-defined lattice BRST 
formulation on the projected manifold and it can serve as a non-perturbative definition of BRST 
symmetry in the continuum limit. 

The Gribov copies of compact U(l) are of course lattice artifacts. However, the gauge group is 
a direct product of odd-dimensional and compact manifolds (circles in this case), as it is for every 
SU(A0 gauge theory on the lattice also. This is sufficient to conclude that it also shares the same 
0/0 problem with S\J(N) theories. Moreover, it turns out that the Neuberger 0/0 problem in SU(A0 
is avoided when that of its maximal Abelian subgroup [/(l)^ -1 is [B], because the coset mani- 
fold consist of even dimensional spheres whose Euler characteristic is 2. It is therefore obviously 
important to understand the Gribov copy and Neuberger problems in compact U(l) gauge models. 

At the same time, for the trivial orbit the functional F in Eq. ( |l.l| ) also represents the Hamil- 
tonian of an XY model which stands alongside the Ising and Heisenberg models as one of the 
most intensely studied systems in statistical mechanics and condensed matter theory. For a random 
orbit, it is that of a random phase XY model (RPXYM). These have been widely used as simple 
models of classical superconductivity, to study aspects of high-r c superconductors, to describe the 
XY magnet with random Dzyaloshinski-Moriya interactions or a positional disordered Josephson 
junction array, to name a few. In either case, the calculations (related to Domain Wall Renormal- 
ization Group studies) boil down to obtaining the global minimum of the Hamiltonian with respect 
to the 0-variables for a given set of 0^'s (see, e.g., [24]). In Ref. [25], it was shown in lower di- 
mensional models that one can first perform a duality transformation then followed by a numerical 
minimization so that a minimization algorithm has to search only within the space of local min- 
ima. Recently, Akino and Kosterlitz [|4|] have used this method to obtain the global minimum of 
a RPXYM Hamiltonian. This method, while giving more confidence in the final results with less 
numerical effort, still relies on numerical algorithms which are known to fail for larger systems. 

The knowledge of all stationary points is important in recent studies of the potential energy 
landscape of the classical statistical mechanical systems such as the RPXYM and its relation to 
phase transitions [20, 21, 22]. An effort to numerically classify all Gribov copies of compact U(l) 
in 2 dimensions was done in [|2"6|]. 

Here, we take up the general problem of finding all solutions to the standard LLG equations 
for compact U(l) or the stationary points of the corresponding XY model. The corresponding 
equations are non-linear which makes them difficult to deal with analytically or numerically. We 
observe, however, that this non-linearity is polynomial-like and the equations can be transformed 
into systems of polynomial equations. Then we can use Algebraic Geometry techniques to solve 
these systems and obtain all stationary points of the corresponding function and in particular the 
global minimum accurately. This procedure is applicable to a wide range of interesting physical 
problems such as obtaining all stationary points of the classical Hamiltonians of Ising and Heisen- 
berg models; solving classical field equations for pure lattice gauge field theories, Abelian Higgs 
model, etc.; finding a generic form of SU(/V) matrices; obtaining the vacua of Supersymmetric po- 
tentials; and many more (see Ref. [jl|]). The corresponding gauge-fixing equations for the modified 
LLG, proposed in Refs. [^, [3]], for the compact U(l) case can also be shown to have polynomial- 
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like non-linearity but their structure is more complicated [JTJ] _ Here, we restrict the discussion to the 
standard LLC 

In the one-dimensional lattice case, the LLG equations have been exactly solved for both peri- 
odic and anti-periodic boundary conditions elsewhere [jl|, ||, ||] where the corresponding equations 
are treated as systems of linear equations. However, for higher dimensional lattices the corre- 
sponding equations are genuinely non-linear and the same method can not be used. We will show 
that the systems of these non-linear equations can be viewed as systems of multivariate polyno- 
mial equations and therefore can be described by the language of Algebraic Geometry. We will 
then discuss two promising approaches to solve the polynomial systems: the numerical polynomial 
homotopy continuation (NPHC) method and the computational Algebraic Geometry approach (in 
the Appendix) and present our results. Though our goal is to solve the gauge fixing conditions for 
higher dimensional lattices, we will explain the methods using the one-dimensional case because: 
(1) we already know the exact solutions for this case and so the results from the new approaches 
will have a precise comparison, (2) the one-dimensional case, interpreted as an Algebraic Geometry 
problem, provides all the essence of the higher dimensional generalizations. 



2. The Extremizing Equations as Polynomial Equations 

Here, we show that the problem of solving the extremizing equations in terms of the 6- 
variables can be transformed into that of solving a system of multivariate polynomial equations. 
We first note that, for a one-dimensional lattice having n lattice sites, the directional index /I is 
irrelevant and ft is just 1. Thus, the corresponding extremizing equations are 

/;-(0 e ) = sin(0 ; e )-sin(0 ; e 1 )=O, (2.1) 
where = (0,- + — 0,) £ {—K, 7t], for all i = 1, .. . ,n. These are the gauge-fixing equations 



(1.2) on the lattice and reproduce the continuum Landau gauge (d^A^ = 0) in the naive continuum 



limit. These equations are also the steady state equations of a variant of the famous Kuramoto 



model to study synchronization in mathematical biology Q23Q. It can be shown that by using anti- 
periodic boundary conditions on both 0- and 6 -variables the global gauge freedom can be com- 
pletely fixed [j^, H|, |l]]. While dealing with periodic boundary conditions one can get rid of the 
residual gauge freedom by taking, say, 6 n to be zero and removing the equation /„ = from the 
system [JIJ]. For simplicity, we take the trivial orbit case (i.e., all 0, = 0). Using the trigonometric 
identity, sin(x + y) = sinjccos y + sinycos* and writing cos 0, = c,- and sin 0, = s,, we get 

fi(c,s) = Ci(s i+ i+Si^i)-Si(ci + ]_ +q_i) = , (2.2) 

for all i = 1, . . . ,n. This is merely a change of notation. However, we add additional equations in 
the system for each site i, namely, 

gi (c, S )=sf + cj -1 = 0, (2.3) 

for all i = 1, ... ,n. Now, the combined system of all fi(c,s) and gi(c,s) is not just a change of 
notation but all the c, and Sj are algebraic variables and the equations are multivariate polynomial 
equations, i.e., the fact that c, and s, are originally sin 8i and cos 0, is taken care of by the constraint 
equations (|Z3|). In general, for the one-dimensional lattice with n lattice sites, we have in total 2n 
polynomial equations and In variables. 
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3. Numerical Polynomial Homotopy Continuation Method 

In the Appendix, a method to exactly solve systems of multivariate equations, called the 
Groebner basis technique, and its current status are briefly explained. Here we discuss a numerical 
method to solve a system of multivariate equations, called the numerical polynomial homotopy 
continuation (NPHC) method. To explain the NPHC method let us consider a system of multivari- 
ate polynomial equations, say P(x) = 0, where P(x) = (pi(x), . . . ,p m (x)) T and x = (x\ , . . . ,x m ) T , 
that is known to have isolated solutions, e.g., the above mentioned LLG equations after eliminating 
the global gauge freedom. Now, there is a classical result, called the Classical Bezout Theorem, that 
asserts that for a system of m polynomial equations in m variables, for generic values of coefficients, 
the maximum number of solutions in C m is fl^Li 4> where d[ is the degree of the ith polynomial. 
This bound, called the classical Bezout bound (CBB), is exact for generic values (i.e., roughly 
speaking, non-zero random values) of coefficients, e.g., for the one-dimensional LLG equations 
with n number of lattice sites and with periodic boundary conditions, this number is 2 2n (because 
there are 2n polynomials each of which is a degree 2 polynomial). The genericity is well-defined 



and the interested reader is referred to Ref. [27, Eq] for details. 



Using the so-called Bernstein-Khovanskii-Kushnirenko theorem, a tighter bound, which takes 
the sparsity of the system into account, can be computed. This bound is called the BKK root count 



or the mixed volume (or stable mixed volume (SMV) Q28| , p2| , [33[ ] in general). The reader is referred 
to the above references to get a precise definition of the mixed volume and its computation. 

Based on either of these bounds on the number of complex solutions, a homotopy can be 
constructed as 

H(x,t) = y(l-t)Q{x)+tP{x), (3.1) 

where /is a random complex number. Q{x) = {q\ (jc), . . . ,q m (x)) T is a system of polynomial equa- 
tions with the following properties: (1) the solutions of Q(x) = H(x,0) =0 are known or can be 
easily obtained. Q(x) is called the start system and the solutions are called the start solutions, (2) 
the number of solutions of Q(x) =H(x,0)=0 is equal to the CBB or SMV for P(x) = 0, (3) the so- 
lution set of H(x,t) = for < t < 1 consists of a finite number of smooth paths, called homotopy 
paths, each parameterized by t E [0, 1), and (4) every isolated solution of H(x,l) = P(x) = can 
be reached by some path originating at a solution of H(x,0) = Q(x) = 0. One can then track all the 
paths corresponding to each solution of Q(x) = from t = to t = 1 and reach P{x) = = H(x, 1). 
Here, a randomly chosen complex y ensures the paths are well-behaved. By implementing an effi- 
cient path tracker algorithm, all isolated solutions of a system of multivariate polynomials system 
can be obtained. 



To illustrate how the method works, we take a univariate polynomial from Ref. [27] to make 
this discussion self-consistent, say z 2 — 5 = 0, pretending that we do not know its solutions, i.e., 
z = izy/S. We first define a family of problems as 

H(z,t) = {l-t){z 2 - l) + t{z 2 - 5) = z 2 - (1 +40 = 0, (3.2) 

where t £ [0, 1] is a parameter. For t = 0, we have z 2 — 1 = and at t = 1 we recover our original 
problem. The problem of getting all solutions of the original problem now reduces to tracking 
solutions of H(z,t) = from t = where we know the solutions, i.e., z = ±1, to t = 1. The choice 
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of the piece z 2 — 1 in Eq. (3.2), called the start system, should be clear now: this system has the 
same number of solutions as the CBB of the original problem and is easy to solve. Now, one of the 
ways to track the paths is to solve the differential equation that is satisfied along all solution paths, 
say z* (t) for the ith solution path: 

dH( z *(t),t) _ dH(z*(t),t)dzttt) , dH(z*(t),t) _ Q _ (3 3) 

dt dz dt dt 



This equation is called the Davidenko differential equation. Inserting (3.2) into this equation, we 
have 

dz*(t) 2 

1 v ' = (3 4) 

dt z*(t) K ' 

This initial value problem can be solved numerically (again, pretending that an exact solution is not 

known) with the initial conditions as z \ (0) = 1 and z\ (0) = — 1 . The other approach is to use Euler 's 

predictor and Newton's corrector methods. We do not intend to discuss the actual path tracker 

algorithm used in practice, but it is important to mention that in the path tracker algorithms used 

in practice, almost all apparent difficulties have been resolved, such as tracking singular solutions, 

multiple roots, solutions at infinity, etc. For the sake of completeness, we should also mention 

here that in the actual path tracker algorithms the homotopy is randomly complexified to avoid 

singularities by y, called the gamma trick, i.e., taking 

H{z,t) = y{\-t){z 2 -\)+t{z 2 -5) =0. (3.5) 



There are several sophisticated computational packages such as PHCpack p9[], PHoM [31 ], Bertini 
[0 and HOM4PS2 [ ^0| ] which can be used to solve systems of multivariate polynomial equations. 
They are all available as freewares from the respective research groups. In the Appendix, we 
compare the Groebner basis technique and the NPHC method with technical remarks. 

3.1 Results 

We present the results by classifying the obtained solutions in terms of the number of positive 
and negative eigenvalues of the corresponding Hessian matrix, or the Faddeev-Popov operator, 
because then we can use the Neuberger zero as a necessary condition for having all the solutions. 

3.1.1 Anti-periodic Boundary Conditions 

We now explore the simplest non-trivial case in higher dimensional lattices which is the stan- 
dard LLG functional on a 3 x 3 lattice with the trivial orbit and anti-periodic boundary conditions 
(i.e., the classical XY model) and also a random orbit case (i.e., randomly chosen ( jU £ (—71,71]). 
For both cases, the CBB was 262144 and the SMV was 148480. For the trivial orbit case, the total 
number of real and complex solutions was 10738 out of which there were 2968 real solutions, i.e., 
Gribov copies. Similarly, for the random orbit case, the total number of real and complex solutions 
was 20558 out of which there were 2480 real solutions 2 . It is important at this stage to mention a 
few specifics about these solutions. We give details for the trivial orbit case, which are similar to 
the other cases, below. 



2 Note that we do not intend to compare the efficiency of the available packages. However, to give an idea how 
computationally intensive such a calculation is, we note that the HOM4PS2 package took around 65 and 120 minutes, 
respectively, to run these systems on a Linux single-processor desktop machine. 
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boundary conditions 


orbit 


real solutions 


sum of 






total number 


zero det. 


non-zero det. 


signs of det. 


anti-periodic 


trivial 


2968 


1152 


1816 







random 


2480 





2480 





periodic 


trivial 


270 


52 


218 







random 


224 





224 






Table 1: An overview on the number of real solutions for the different cases 
(orbits and boundary conditions). See Section 3 for the detailed description. 



Case 


i 





1 


2 


3 


4 


5 


6 


7 


8 


9 


antiperiodic b.c, trivial 


Ki 


2 


18 


216 


342 


330 


330 


342 


216 


18 


2 


antiperiodic b.c, random 


Ki 


2 


58 


202 


402 


576 


576 


402 


202 


58 


2 


periodic b.c, trivial 


Ki 


8 


28 


38 


42 


56 


38 


7 


1 






periodic b.c, random 


Ki 


2 


12 


30 


61 


72 


38 


8 


1 







Table 2: Summary of the number of real solutions Ki with i negative eigenvalues 
on the 3x3 lattice, for the different cases (orbits and boundary conditions). 



1. A solution means that a set of values of Sj's and c,'s satisfies each of the 18 equations with 
tolerances 1 x 10~ 10 . All the solutions come with real and imaginary parts. A solution is a 
real solution if the imaginary part of each of the 18 variables (i.e., all Sj and c;) is less than or 
equal to the tolerance 1 x 10~ 6 (below which the number of real solutions does not change, 
i.e., it is robust for all the cases we consider in this discussion). The original trigonometric 
equations are satisfied with tolerance 1 x 10~ 10 after s, and c, are transformed back to 0,. 
All these solutions can be further refined with an arbitrary precision. This is a remarkable 
success of the method because then these solutions are close to the exact solutions. 

2. The sum of the signs of the Faddeev-Popov determinants for all real solutions for the trivial 
and random orbit cases is ~ 10~ n , which is numerically zero, yielding the expected Neu- 
berger zero as mentioned in the Introduction. 

3. For the trivial orbit case, there are exactly 1 152 real solutions which have zero Faddeev- 
Popov determinant with tolerance 1 x 10~ 8 . These solutions constitute the set of singular 
loci or the Gribov horizons. They can be further classified in terms of the number of zero 
eigenvalues of the Faddeev-Popov operator at each of the solutions. This amounts to clas- 
sifying the singular solutions of the polynomial system in terms of their multiplicities using 
the so-called deflation singularities technique^]. For both the random orbit cases, there is 
no Gribov horizon, i.e., all the solutions are non-singular. 

4. For the trivial orbit case, the remaining 1816 real solutions have nonsingular Faddeev-Popov 
determinant, i.e., nonsingular solutions. All the nonsingular solutions of both the trivial and 
the random orbit cases can be classified by the number of negative eigenvalues (Table ^|), 
which shows the expected two-fold symmetry giving rise to the Neuberger zero. 

See Table for the summary. 
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3.1.2 Periodic Boundary Conditions 

After eliminating the constant zero mode for the 3x3 lattice case with periodic boundary 
conditions (by taking, say, two of the 0,'s at the corners of the lattice to be 0, without losing 
generality See Ref. [jl|] for more details), we have 14 equations and 14 variables. The number 
of solutions is 270 out of which 218 are non-singular and 52 are singular. The Faddeev-Popov 
operator is now a 7 x 7 dimensional matrix. The sum of signs of the Faddeev-Popov determinants 
is zero. Similarly, for a random orbit the number of real solutions is 224 and again the sum of signs 
of the Faddeev-Popov determinants is 0. The results are summarized in Table |2| 



4. Summary 

By taking the standard lattice Landau gauge-fixing (LLG) functional (the classical Hamilto- 
nian of the random phase XY model (RPXYM)) as an example, we showed that the non-linearity 
of the corresponding gauge-fixing equations (extremizing equations) is polynomial-like. With ad- 
ditional constraint equations, the combined system of equations can be treated as a system of poly- 
nomial equations with all variables defined over C. Though we are interested in only real solutions 
of these polynomial systems, we used Complex Algebraic Geometry concepts over Real Alge- 
braic Geometry due to the stronger results available in the former. To solve the corresponding 
polynomial systems exactly, one can use an elegant algorithm, called the Buchberger algorithm 
(see Appendix A). Though the algorithm works well for the one-dimensional lattice, for the 3x3 
lattice, it is beyond the capabilities of a standard desktop machine. However, the method can be 



useful in finding a parametrization of the SU(/V) matrices [38] and hence can then be used, for 
example, to extend the method of restricting a path integral of the SU(A0 Georgi-Glashow model 
to the 't Hooft-Polyakov sectors by using the twisted C-periodic boundary conditions [Q] for odd 
N and to precisely define the modified lattice Landau gauge for the SU(3) case. 

For the two-dimensional case, we used the Numerical Polynomial Homotopy Continuation 
(NPHC) method, which gives all solutions of a polynomial system numerically. This method 
does not suffer from the technical difficulties of the Buchberger algorithm or its variants, and in 
principle this method can find the solutions for those systems which may be intractable for the 
former method. By computing the eigenvalues of the Hessian matrix at these solutions, one can 
obtain the global minimum, up to machine precision, of a multivariate function whose extremizing 
equations can be translated to polynomial equations. 

The 3x3 anti-periodic boundary conditions system for the trivial and a random orbit can be 
solved using the NPHC method and that there is an exact cancelation of the signs of the Faddeev- 
Popov determinants in this case, i.e., the Neuberger zero. We classified all the solutions in terms 
of the number of negative and positive eigenvalues of the corresponding Faddeev-Popov operator. 
Apparently, the number of Gribov copies is orbit dependent, though the Neuberger zero is orbit 
independent. It would be important to compute the number of solutions for different random or- 
bits for the 3x3 lattice with the modified LLG using the NPHC method. Efforts to solving the 
corresponding gauge-fixing equations for the modified LLG for compact U(l) and also the linear 
covariant gauge-fixing equations for compact U(l), in the spirit of [4-0, jUl], are in progress. It 
would also be very interesting to solve the Langevin dynamics equations for the XY model and 



other models [39] using algebraic geometry methods. 
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A. Computational Algebraic Geometry 

Since many important and useful results in algebraic geometry are mainly developed for the 
complex field C, we treat our polynomial equations as having complex variables, and in the end 
we use only real solutions and throw the remaining solutions away (see, e.g., Jl8|]). One can then 
use an important result from computational algebraic geometry: after specifying an ordering of the 
monomials of the system (denoted as x >- y >~ z for the ordering in which x is placed prior to y 
and y is placed prior to z), one can transform a given system of multivariate polynomial equations, 
referred to as an ideal, to another one which has the same solutions but is easier to solve. This new 
system is called a Groebner basis for the given monomial ordering. 

There is a well defined procedure to find a Groebner Basis for any given ideal and monomial 
order, called the Buchberger algorithm. It should be noted that the Buchberger algorithm reduces 
to Gaussian elimination in the case of linear equations, i.e., it is a generalization of the latter. Simi- 
larly it is a generalization of the Euclidean algorithm for the computation of the Greatest Common 
Divisors of a univariate polynomial. Recently, more efficient variants of the Buchberger algorithm 



have been developed to obtain a Groebner basis, e.g., F4, F5 and Involution Algorithms [34]. Sym- 
bolic computation packages such as Mathematica, Maple, Reduce, etc., have built-in commands to 
calculate a Groebner basis for a given monomial. Singular, COCOA and McCauley2 are special- 
ized packages for Groebner basis and Computational Algebraic Geometry, available as freeware. 
MAGMA is also such a specialized package available as a non-free package. Rather than going 
into the details of the specifics and technicalities of this algorithm, we dive into the practical appli- 
cations of the Groebner basis technique relevant to our problem and refer the reader to the above 
mentioned references for further details. 

A.l The Groebner Basis Technique At Work 

Here, we provide a practical example of how the Groebner basis technique can be used for 
our systems. Firstly, for the polynomials for the one-dimensional lattice with three lattice sites, 
anti-periodic boundary conditions and the trivial orbit case, the corresponding ideal (denoted as the 
polynomials put between '<' and '>') is 

h = < -C2S1 -C 3 Si+CiS 2 -CiS 3 ,C 2 Sl -CiS 2 -C3S2+C2S3,-C3Sl+C 3 S2-CiS3-C2S3, 

c? + 4-l,cl + ji-l,cf + j|-l>. (A.l) 
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For the ideal in Eq. ( A. 1 ), a Groebner basis for the so-called lexicographic ordering c\ >~ s\ >- C2 >- 

$2 >- ^3 y S3 is 

G3 = < -S 3 +sl,C3S3,-\+cj+sj,-S2 + S 2 sj,C3S2,sl-sj,C2S3,C2S2, 

-1 +c\ + sl,-Si +Sisj,C3Si,C 2 Susf -sl,C\S3,ClS2,ClSl,-l +c\ + sj > . (A.2) 

As noted earlier, the solutions of this system are the same as the original system. Here, the first 
equation in G3 is a univariate polynomial in variable S3 and solving it is simple because it can be 
factorized as S3 (^ — 1 ) =0 giving S3 = and S3 = ±1. Using this and back-substitution, all the 
solutions can be obtained and are as follows: 

(cuawus^) = {(0,0,0,±1,±1,±1),(±1,±1,±1, 0,0,0)}. (A.3) 

Thus, the solution space of the system G3, called the variety of G3 and denoted as V(G^), is the 
above mentioned set of 16 isolated points in a 6-dimensional affine space. This is the known result 
from Refs. ^[ [I]] in terms of 0,'s, i.e., 

(01,02,03) = {all 2 3 permutations of 0/7l}U {all 2 3 permutations of -f/f }. (A.4) 



All of these solutions are correctly reproduced in Eq. (A.3). Similar computation can be performed 
for the periodic boundary conditions case 

For an ideal which is known to have only isolated solutions (called a zero-dimensional ideal) 
with a lexicographic ordering one can always find a Groebner basis in an upper diagonal form, 
analogous to the Gaussian Elimination method, such that at least one polynomial is univariate with 
the others having an increasing number of variables. 

Interpreting the gauge fixing equations in terms of Algebraic Geometry allowed us to deal 
with the actual non-linearity of the equations, rather than treating them as linear equations as done 
in Refs. [pi 0, [3|]. Specifically, in this interpretation the corresponding method does not make any 
distinction between the equations arising from a one-dimensional lattice or those arising from a 
higher dimensional lattice. In theory, as long as one can obtain a Groebner basis, the equations can 
be exactly solved. For the one-dimensional case, it was easy to go up to a fairly big number of 
lattice sites on a single machine. However, in general, obtaining a Groebner basis is very difficult 
due to an algorithmic complexity, known as Exponential Space complexity, which roughly means 
that the RAM required by the computation blows up exponentially. In particular, on a regular single 
desktop machine with 2 GB RAM, we could not obtain a Groebner basis for the trivial orbit case 
(i.e., the classical XY model) and with anti-periodic boundary conditions for a 3 x 3 lattice, using 
Singular3.2. The corresponding system is made of 18 equations, each of degree 2, in 18 variables. 
However, a more powerful machine should be able to obtain a Groebner basis for these systems. 

Recently, V. Gerdt and Daniel Robertz were able to compute a Groebner Basis for this system 
over F2 = {0,1} using a Linux machine with AMD Opteron (TM) processor 285 2600 MHz, 4 
cores, 16 GB memory and with Magma (V2. 14.14) in less than 10 hours with degree reverse 



lexicographic ordering [35]. This is a very important step towards solving the system. However, 
since the computation is over F2 not all the solutions over C 18 can be obtained. 
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A.2 Comparison between the Groebner basis technique and the NPHC Method 

The NPHC method is strikingly different from the Groebner basis technique in that the algo- 
rithm for the former suffers from no known major complexities. Moreover, the path tracking is 
embarrassingly parallelizable, because all the start solutions can be tracked completely indepen- 
dently of each other. This feature along with the rapid progress towards the improvements of the 
algorithms makes the NPHC well suited for many physical problems arising in condensed matter 
theory, lattice QCD, etc. 

For most systems of polynomial equations in practice, we do not know the actual number of 
solutions from the beginning. So even after obtaining all solutions through the homotopy continua- 
tion method, we would still prefer to have some kind of verification of the solutions. In the standard 
LLG fixing case, we used the Neuberger zero as a necessary condition. However, this is certainly 
not a sufficient condition and further checks may be required for bigger systems. This is called 
certifying the solutions, which has been recently developed using the so-called Tropical Algebraic 



Geometry ||37|]. We anticipate that our results mentioned above will serve as ideal test systems for 
all such new ideas. 

To solve the corresponding equations for the periodic boundary conditions case, another method 
is to leave the constant zero modes in the system and use the Numerical Algebraic Geometry 
(NAG) p2|], which is the generalization of the NPHC method to positive-dimensional ideals. 

Recently, a computational Algebraic Geometry approach to find vacuum configurations in 
string phenomenology (where mathematically the problem reduces to minimization of the so-called 



superpotentials which exhibit polynomial-like non-linearity) is proposed in Refs. [ jig , |Psfl . A very 
efficient Mathematica package, called STRINGVACUA, has been developed with an interface with 
the computer algebra system Singular. This package relies on the Groebner basis technique in solv- 
ing the extremizing polynomial equations. It is anticipated that incorporating the NPHC method 
can extend the applicability of the package and the approach beyond its present status. 

It should be noted that there are several promising methods such as the Discriminant Vari- 



ety [|43|] and finding real solutions out of complex curves Q44Q from computational real algebraic 
geometry which are emerging as alternatives to the above mentioned methods in some of the im- 
portant problems mentioned in the Introduction. 
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